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We study the asymptotic position distribution of general quantum walks on a lattice, including 
walks with a random coin, which is chosen from step to step by a general Markov chain. In the 
unitary (i.e., non-random) case, we allow any unitary operator, which commutes with translations, 
and couples only sites at a finite distance from each other. For example, a single step of the 
walk could be composed of any finite succession of different shift and coin operations in the usual 
sense, with any lattice dimension and coin dimension. We find ballistic scaling, and establish a 
direct method for computing the asymptotic distribution of position divided by time, namely as 
the distribution of the discrete time analog of the group velocity. In the random case, we let a 
Markov chain (control process) pick in each step one of finitely many unitary walks, in the sense 
described above. In ballistic order we find a non-random drift, which depends only on the mean 
of the control process and not on the initial state. In diffusive scaling the limiting distribution 
is asymptotically Gaussian, with a covariance matrix (diffusion matrix) depending on momentum. 
The diffusion matrix depends not only on the mean but also on the transition rates of the control 
process. In the non-random limit, i.e., when the coins chosen are all very close, or the transition 
rates of the control process are small, leading to long intervals of ballistic evolution, the diffusion 
matrix diverges. Our method is based on spatial Fourier transforms, and the first and second order 
perturbation theory of the eigenvalue 1 of the transition operator for each value of the momentum. 
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I. INTRODUCTION 



A quantum walk, the counterpart of a classical random walk, models a quantum particle moving randomly in 
discrete time steps on a lattice, but contrary to the classical case one has to consider a quantum particle with an 
internal degree of freedom, usually called the coin space, in order to get nondeterministic behavior [TJ [5]. Frequently, 
the dynamics of a quantum walk is decomposed into a unitary operation acting on the internal degree of freedom, 
called the coin operation, followed by a state dependent spatial shift operation. This constructive definition 
permits many interesting examples of quantum walks and moreover it yields a decomposition into experimentally 
realizable operations, which already have been implemented [5]. However, in general there is no need to restrict 
to dynamics that can be decomposed into coin operations and conditional shifts. Indeed, we will consider a more 
general axiomatic definition of quantum walks without any intrinsic connection between the lattice dimension and 
the number of internal states of the particle, which is suggested by the constructive approach to quantum walks. 
Without adhering to such a connection, quantum walks with memory as introduced by McGettrick [1] can easily be 
formulated in our language and their analysis, in particular the calculation of their asymptotic behavior, is much 
simplified. The definition of quantum walks we are going to incorporate imposes two axioms on the time evolution. 
The first one is translation invariance, i.e. the time evolution is homogeneous in space, meaning it commutes with 
lattice translations. The second axiom is a locality condition, we want to exclude infinite propagation speed of the 
particle and therefore we assume that the maximal step size of the particle is bounded. Interestingly, it turns out 
that for one dimensional lattices and unitary time evolution this axiomatic definition exactly matches the constructive 
approach J5], whereas in higher dimensions this is not clear. 

In this paper we consider quantum walks on lattices in any space dimension. Their internal degree of freedom is 
described in any finite dimensional Hilbert space. According to our axiomatic definition the dynamical step will be 
translation invariant and the step size can be many lattice constants long but finite. Another assumption we may 
impose on the time evolution is unitarity. Although we also discuss general decoherent quantum walks compatible with 
our definition, we will later consider a particular class of decoherent quantum walks emerging from a set of unitary 
ones. What we want to study is disorder in time, rather than in space: We assume that there is some collection of 
possible steps, from which a classical Markov process chooses one instance. Thus successive steps are not independent, 
although steps separated by a long time will be practically independent. 

For this class of models we aim to compute the asymptotic probability distribution for the position variable Q(t) 
after t steps. Such analysis has been carried out before, but to the best of our knowledge not in full generality. Most 
of the analysis so far is concerned with unitary quantum walks [TJ O [7] or under certain constraints, e.g. in one 
lattice dimension and with step size one [SlllO] . A numerical study of certain examples of unitary quantum walks in 
two lattice dimensions was performed by Mackay et al. [H] and Kollar et al. [12]. Recently Baryshnikov et al. [13] 
and Bressler et al. I14j computed the limiting probability distributions for various one and two dimensional unitary 
quantum walks. Also, the study of quantum walks subjected to decoherence is mostly concerned with particular models 
or subclasses of quantum walks, see e.g. the review article by Kendon |15j and references therein for an overview on 
decoherent quantum walks on one dimensional lattices, cycles or hypercubes. The kinds of decoherence considered 
can be roughly classified as measurement induced decoherence [TBH^D] and decoherence caused by sloppy control of 
quantum operations |21H24j . A subclass of the second decoherence mechanism which we are going to consider is the 
case when in each time step a unitary quantum walk is chosen from a given set with a certain probability, for example 
the coin operation of a one dimensional quantum walk may be chosen from a fixed set of unitaries at each time step 
[25H28] . A more general kind of decoherence compatible with our definition of a quantum walk was introduced by 
Aimabestani et al. [29] and Annabestani et al. |30] , and an analysis of the asymptotic behavior was performed based 
on the calculation of the first and second moment of the distribution. 

One key question we address is whether the spreading is ballistic, i.e. whether Q(t)/t converges in distribution, 
or whether it is diffusive, i.e. Q(t)/y/i has a meaningful limit as t — » oo. Now without randomness, i.e. when we 



use always the same step, it is known that the spreading is ballistic. In this case, as we show in Sec. IV A the 
asymptotic distribution of the particles position is determined by the dispersion relation of the unitary quantum 
walk, more precisely, its derivative with respect to momentum, the group velocity, dictates the asymptotic behavior. 
Commonly, the asymptotic distribution shows characteristic peaks, which can be understood as caustics of the 
dispersion relation, see Fig. [I] Such caustics have also been observed by Baryshnikov et al. [13 , we give a detailed 
description of this concept in Sec. |IV A| Since a classical random walk satisfies our description it is also clear that 
for stringent randomness we expect diffusive scaling. In fact, diffusive behavior is a common feature of decoherent 
quantum walks, e.g. for quantum walks with multiple coins [2"2l |2"5H2"T] or dynamic gaps in the lattice [T71H2GIT] this 
has been observed. On the basis of computations of the second moments, it has been shown recently [29] that diffusive 
scaling holds even for very low randomness, i.e. when either all the steps used are nearly the same, or if one step is 
chosen almost always. Our results support this conclusion. In addition, we compute the asymptotic distribution of 
Q(t)M for every initial state. It is Gaussian in every momentum component (in a sense specified below). As is to 



FIG. 1. The figure shows a caustic of a unitary quantum walk on a two-dimensional lattice. Closed curves correspond to the 
image of the map p h-> Vcj(p) applied to a discrete set of coordinate lines in momentum space [0,27r) 2 . The caustic, i.e. the 
region where the line density diverges, is exactly the region where the asymptotic probability distribution of the quantum walk 
exhibits peaks. 



be expected, the diffusion constant (the limit of (Q 2 (t))/t) diverges at low randomness. In such cases the system will 
initially evolve ballistically, and then exhibit a crossover to diffusive scaling, see Fig. [T] and [2j 

The methods used so far to obtain the asymptotic position distribution can be summarized as Fourier methods, 
see e.g. Ambainis et al. PQ, Grimmett et al. [7j, and combinatorial calculations as used by Konno Another 
approach, closely related to Fourier methods is the generating function formalism |10[ I19j . Our method involves 
Fourier methods as a way of describing a translation invariant system and its time evolution. The second main 
ingredient of our formalism is perturbation theory. We phrase the problem of finding the asymptotic distribution of 
the scaled random variables Q(t)/t and Q(t)/\/i as the problem of applying a perturbed operator infinitely many 
times to an eigenvector of the unperturbed operator in the limit of vanishing perturbation. This becomes a nontrivial 
problem because both limits do not commute. The advantage of our method is that it gives us full information about 
the limiting position distribution of a large class of quantum walks, whereas other approaches often only give partial 
information, e.g. about the first two moments of the limiting position distribution. 

Our paper is organized as follows: In the next section we will introduce notation and basic concepts. We close 
by describing the most general translation invariant quantum evolution with strictly finite step sizes, extending a 
result prematurely claimed to be exhaustive by Annabestani et al. [29] . The type we actually consider in most of the 
paper is a special case and essentially the same as the one considered by Annabestani et al. [29] ■ We then briefly 
review the non-random case, which is a prerequisite to understanding the general case. In Sect. V A VB| and |VD| 
we then state our result with all assumptions spelled out, and with a description of the procedure to compute the 
limit distributions. In the final section we gather some examples of quantum walks which violate our assumptions 
and analyze their asymptotic behavior. 



II. THE SYSTEMS 



The underlying lattice of the system will be denoted by X, and is a subset of K s . We will always take X = Z s , 
but we would like to stress that this covers also much more complex periodic structures, like triangular lattices, 
honeycombs, Kagome lattices, and so on in higher dimensions. Indeed, for us the role of the lattice is mainly that of 
the abstract symmetry group. For a general periodic pattern, the translation symmetries are by vectors of the form 
Si=i X A f° r some basis of lattice vectors and integer coefficients Xj. One can then choose a unit cell, so that 
every point in M. s is uniquely obtained as a point in that cell plus a lattice vector. Then we consider all lattice points 
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(c) (d) 

FIG. 2. (color online) Figures (a)-(c) show the time evolution of a quantum walk with index three and three internal states 
for 50, 150 and 300 time steps. The green/dot-dashed curve corresponds to the undisturbed walk with a coin that is constant 
in time, whereas the red/solid and blue/dashed curves correspond to a Bernoulli process where in every time step instead of 
the coin of the green/dot-dashed quantum walk a different coin is applied with probability 0.01 or 0.2, respectively. Graph (d) 
shows the standard deviation of the position distributions with a visible crossover from ballistic to diffusive spreading behavior 
for the perturbed quantum walks. 

inside the unit cell as one super-site with internal structure, and take the integer coefficients xi themselves to label 
the lattice translations. This brings us back to the choice X = Z s . The position degree of freedom of the walking 
quantum particle is thus described in the Hilbert space t 2 (X), the square-summable functions ip : X — >• C. 

The internal states of the walking particle are described by a finite dimensional Hilbert space JC. In many papers 
on quantum walks this is called the coin space. We will follow this terminology although the dynamics we consider 
usually does not have a simple decomposition into coin operation and conditional shifts. The Hilbert space of the 
system is thus i 2 (X) g) tC, and it will be convenient to identify this with £ 2 (X, JC), the set of functions ip : X — > JC 
such that H^ll 2 = J2xex \\*P( X )\\ 2 < 00 • Translations now act as shifts of the argument, i.e. (U x ip)(y) — ip{y — x). 

Since we look at translation invariant systems, the analysis can be simplified considerably by taking Fourier trans- 
forms. Let X denote the dual group, in our case concretely parameterized as the space of momentum vectors 
p £ (— 7r, ir] s C M. s . The Fourier transform is then the unitary operator 

(^)(p)-(2^)- s / 2 ^e^(x). (1) 

Thus J 7 is a map from £ 2 (X,JC) to £ 2 (X,JC). We can also think of this equation as the representation of a gen- 
eral function on M. s , which is periodic with respect to the shifts in (2irZ) s . Of particular importance will be the 
trigonometric polynomials given by such a series with just finitely many non-zero terms. 

A unitary quantum walk is given by a unitary operator W, which commutes with translations and has the 
property that (Wip)(x) depends only on the values of ip(y) such that x — y £ N , where Af is a fixed finite set called 
the neighborhood scheme of W. Because W commutes with translations it can be diagonalized jointly with the 
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translation operators, i.e. it becomes a multiplication operator in momentum space: 



W)(p) = W r (p)W)(p). (2) 

where, for each p, W(p) is a unitary operator on /C, and each matrix element of W(p) is according to the locality 
condition a trigonometric polynomial. More precisely, the only monomials e tx ' p appearing with non-zero coefficients 
are those with x £ J\f. In electrical engineering such operators are called para-unitaries. One can show [33 35 that, 
at least for s = 1, every such operator can be decomposed into a finite product of p-independent unitaries ("coin 
tosses") and diagonal unitaries which only have elements of the form e lx ' p on the diagonal ("conditional shifts"). 

In order to characterize unitary quantum walks it is useful to consider the determinant of the walk operator W(p), 
which is also a trigonometric polynomial. Since 1/ det W(p) — det(W(p)*) is also a polynomial, it must be a monomial 
e ix-p f or some x £ X. We call this x the index of W and by definition we have 

det W(p) = det W(0) • e ind w ' p . (3) 



Turning now to irreversible processes, let us first introduce the kind of randomness, which is analyzed in our main 
result. We consider unitary walks, such that the particular unitary applied at time t is chosen randomly. In this 
way we want to model experiments in which the randomness comes from fluctuations in some external parameters 
controlling the quantum operations. If the time scale of these fluctuations covers several steps, it is unreasonable to 
assume that the unitaries in successive time steps are independent. Therefore we allow the external parameters to 
be given by a Markov process, which we call the control process. We fix some dependence 7 H> W 1 of the walk 
unitaries on the value 7 of the control process. Then the following two steps are iterated: At time step t the walk 
unitary W 1 {i) acts on the quantum system, where 7(f) is the current value of the control process. Then the next 
value j(t + 1) is drawn according to the transition probability law of the control process. Our aim is to derive the 
long time behavior of this scheme, particularly the distribution of the particle's position at large times. 

Depending on the transition probabilities, we can describe systems in which the walk unitaries of successive steps are 
either strongly correlated, or nearly independent. We assume that in the long run the process goes to an equilibrium 
probability density. This stationary distribution will be taken as the initial distribution. In this way we express 
the condition that there are no correlations between successive runs of the walk, as the statistical data are collected. 
Otherwise, non-trivial correlations would exist between successive runs, and the experimental results would depend 
on how quickly the next run is initiated. Of course, this kind of correlations is possible, and may actually occur 
in experiments, but it would mean a deviation from the paradigm of statistical data collection. Theoretically the 
appropriate response would be to describe not single runs, but batches of runs with controlled relative timing. This 
is a complication we do not want to consider in this paper. So, as already stated, we will assume that the initial 
distribution is the stationary one. 

Of course, further generalizations are possible. Firstly, we may have decoherence in the individual steps, so that 
the single step is not given by a random unitary, but by a completely positive map, whose Kraus operators are not 
multiples of unitaries. We will analyze the general form of such processes in the next section. We will also point 
out (Sect. VAl the asymptotic results which can be obtained for this more general class of decoherent walks using 
the method of the main theorem together with an assumption concerning the general form of the Kraus operators. 
However, in the whole paper we will keep translation invariance. When this is broken, perhaps randomly, a whole new 
range of phenomena appears. The most interesting is the analogue of Anderson localization: Whereas the decoherence 
studied in this paper slows the spreading of the walk from ballistic scaling (~ t) to diffusive scaling (~ i 1 / 2 ), random 
space dependent coins typically stop the spreading altogether (~ t ) in the sense that an initially localized state will 
remain finitely localized for all times with arbitrarily high probability. Here the localization region depends on the 
initial state and the allowed error probability, but not on time. This has been demonstrated by Ahlbrecht et al. [37] 
using methods adapted from the much better studied case of continuous time Anderson localization. Eventually, it 
will be possible to combine randomness in space and time, but this is definitely beyond the scope of the present work. 



A. General form of decoherent, translation invariant quantum walks 



Turning now to decoherent dynamics we will no longer consider pure quantum states ip G £ 2 (X,IC) but density 
operators p, i.e. bounded and positive operators with unit trace. Hence, p is an element of the space B(H) of bounded 
operators on the Hilbert space H — i 2 {X, JC). Generically, a density operator p will be non-translation invariant, e.g. 
p may be supported on finitely many sites of the underlying lattice. Hence, the Fourier transform of p is not given 
by a multiplication operator as in p|, but an operator depending on two variables p and p' . An example of such an 
operator is a state which is a mixture of finitely many pure states ip that are all supported on finitely many sites of 
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the lattice: 

A> = X>|lfcWi| p{p,p')=F P F f = Y,Mi>i(p))^i(p')\ 

i i 

If we consider p(p,p') as a linear map its action and the trace are given by the formulas 

(Tpip)(p) = / dp' p{p,p') ■ ip{p') and trp= / dp tr p(p,p) = dp trp(p) 



where we abbreviated p(p,p) — p{p)- A general non-unitary time evolution of a quantum system is described by the 
concept of a quantum channel W. In the following we will describe the dynamics in the Heisenberg picture, i.e. we 
will evolve the observables of the system rather than the states. Hence, a quantum channel W : B^H.2) — > B(Hi) is 
formally given by a linear map which is completely positive and unital, i.e. for n E N and A E B(H2 <8> C") 

A > => W ® ido(A) > and W(1«J = 1-^ . 

Here, Hi is the Hilbert space of the initial system, and H2 describes the system after the time evolution. The 
representation theorem by Stinespring [38 , 39. states that a completely positive and unital map W can be written 
explicitly as 

W(A) = V*{A® 1 V )V ,VAeB(H 2 ), 

where V : Hi — > %2®T> is an isometry, i.e. V* V = 1-h 1 , and V is called the dilation space. By choosing an orthonormal 
basis e, in V and writing 1-p = J^i l e i)( e il one obtains the Kraus representation [33J 3D] of a quantum channel W 
which reads 



W(A) = K*AK, VA e B{H 2 ) 



The relation between the Ki and the isometry V is given by 

(<t>\Kiil>) ={<}>® e,|V^) ^%,^e%. 

The question we are going to address now is which isometries V, respectively Kraus operators Ki, represent a trans- 
lation invariant quantum walk W. We denote translations by lattice vectors x G X by r x , that is, t x e B(£ 2 (X) ® /C) 
is defined via 

T x {\y){z\®M) = \y + x){z + x\®M , y,z E X , M E B(/C) . 

With this definition translation invariance of the quantum walk W is expressed by 

W(r x {A))=W(A), VxeX,AeB{£ 2 (X)<Z>1C). (4) 

In order to exclude infinite propagation speed we also impose the following locality condition on the quantum walk 
W, which by Q can be chosen translation invariant. We assume there exists a finite neighborhood scheme Af C X 
such that for arbitrary internal states <fi, ip E K, and M E B(JC) 

{k®(j)\W{\x){y\®M)\l®il))=Q if k-x£ Af or y-l £Af. (5) 

The following theorem characterizes all translation invariant quantum walks with finite propagation speed. 

Theorem 1 Let W be a translation invariant quantum walk on % — £ 2 (X) <X> K., that is, W : B(fH) — > B{Ti.) is a 
completely positive map respecting and Then, there exists a dilation space T> and a unitary representation 
{U x }xeX °f X on T> together with operators v y : K, — > K, ® T>, y E Af, such that an isometry V : % — > % ® T> 
representing W is given by 

V\x® 4>) = ^2 \x + y) ® ® U x )v v <j>) , xEX,(j)E]C. (6) 

Conversely, a unitary representation {U X } X £X of X onT> together with operators v y : K, — > K, ® T>, y E Af ', satisfying 
the normalization condition 

v* +y (l K ®& x )v y = l K -6 x0 , (7) 
defines a translation invariant quantum walk W via Q . 
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Corollary 2 By choosing orthonormal bases ej £ D and e a £ K we get the Kraus operators Kj £ B(£ 2 (X) ® JC) 
corresponding to the isometry V of Theorem [7] via 

(z® ep\Kj\x® e a ) = {ep ® e \{l®U x )v z - x \e a ) . 

Proof of Theorem^ 

A general isometry V is given by the relation 

V\x ® 0) = ^ |a; + y) ® Ms) 0) , 

with operators v y (x) : JC — > /C ® 2?. Invariance under translation by i € X requires VE4, =U X ® U X V, where ET^ is 
the translation operator on H and the U x form a representation of X on 2?. The intertwining relation of V leads to 
v v (x) = U x v y (Q) =: U x v y . The locality condition (|5j) assures v y = if y ^ A/" and the normalization condition is a 
consequence of the isometry condition 5 xy ■ (<j)\ip) = {x <g> 0|V*V|y ® V)- ' 

If the spectrum of the operators f7 x is only pure point it is easy to define Fourier transformed versions of V and Kj . 
By choosing an orthonormal basis ej £ V, j £ J, of eigenvectors of the U x and writing U x = Yl eiqj x \ e j)( e j\ with 

qj £lwe get the following corollary. 

Corollary 3 If the spectrum of the U x in Theorem [7] consists only of pure point spectrum and ej labels a common 
eigenbasis of the U x with eigenvalues e lqj ' x , then the isometry V and the Kraus operators Kj in momentum space are 
given by 

(0® ej -|(VV)(p)) = £ e^i^Vy^p + qj)) , cj ) ,yj£lC 

mK^){p)) = E e iM <4KV(p + qj)) , 0, V e £ . 

In general, the operators may also exhibit continuous spectrum. A simple example of such a translation invariant 
quantum walk is the following. Consider TL = ^ 2 (Z), that is, a particle with no internal degree of freedom moving on a 
one dimensional lattice. As a further simplification we will assume M — {0}, i.e. the particle is not moving at all. We 
choose a representation {U x } x ^z of Z on some infinite dimensional dilation space T> and an operator vo : C — > £ 2 (Z). 
The operator V defines a translation invariant quantum walk if the normalization condition ([7]), which reads 

v*(l K ®U x )v =5 x0 , VxeZ, 

is fulfilled. By definition v x is zero if x ^ 0, hence, vq is an isometry. The action of the operator W corresponding to 
V is 

W(\x){y\) = (vo\U y -*\vo) ■ \x){y\. 

Obviously, W satisfies the locality condition ([5| with Af = {0}, it leaves diagonal elements invariant and off-diagonal 
elements are damped exponentially in the number of time steps. A possible choice for U x with continuous spectrum 
is T> = ^(Z), any normalized vector vq £ ^(Z) and the operators U x as shift by x lattice sites. 

General translation invariant quantum walks with momentum transfer are difficult to handle. We must assume 
in the following that there is no momentum transfer. An example where we drop this assumption can be found in 
section IVIDl 

Assumption 1 There is no momentum transfer, i.e. U x — 1, Vx £ 1/ . The isometry V and Kraus operators Kj are 
given by 

(0|(A^)(p)) = £ e^®e>^(p)) , 4>, $ etC, 
where ej labels an orthonormal basis ofT>. 



III. ASYMPTOTIC POSITION BY THE PERTURBATION METHOD 



From the beginning of quantum walk theory the question of asymptotic behavior of the position Q(t) at large 
times t has been one of the main themes. Early papers were to some extent misguided by the analogy with random 
walks and long combinatorial computations of matrix elements were done to evaluate just the special case of a 
Hadamard walk starting from the origin. The physicists in the community quickly brought to bear Fourier methods, 
and these emphasized the analogy not with classical random walks, but with the free particle under a continuous time 
Schrodinger time evolution. In particular, this brought in dispersion relations w(p), and group velocities Vw(p) as the 
relevant quantities in the unitary case. Thus, for a general unitary walk and arbitrary initial state, the computation 
of the asymptotic distribution of Q{t)/t became a straightforward evaluation of expectation values (see below). Initial 
studies on decoherent walks were often limited to a very special noise model, and have almost exclusively considered 
the first and second moments of position. While this is already good enough to distinguish ballistic from diffusive 
transport, it usually remained open how to compute the asymptotic distribution of Q{t)/^fi or, indeed, how to decide 
whether this quantity had a limit distribution. 

We will therefore begin by showing how to focus on the entire distribution of Q from the outset. This will establish 
the perturbation theory of the eigenvalue 1 of the transition operator as the key tool in the further analysis. In 
that introductory section we ignore the control process in order to keep the notation simple. We then look at the 



unitary case (Sect. IV A I, where the first order perturbation theory of a degenerate eigenvalue determines the group 
velocity operator and hence the ballistic scaling. This is then extended to more general processes, including externally 
controlled ones (Sect. [VA ). We then come to our main result, the asymptotic formulas for Markov controlled coined 



walks, in diffusive scaling (Sect. VB ). Finally, we consider the simplified case where the Markov process is of Bernoulli 
type (Sect. VD I, i.e. its transition rates are independent of previous time steps, which means the quantum operations 
are drawn independent and identically distributed in each time step. 

A compact way to characterize a probability distribution of a real vector valued random variable Q £ K s is in terms 
of its characteristic function 

C Q (A) = <e <A -«>. (8) 

Here the bracket denotes expectation, A is a real vector of the same dimension s as Q, and the product in the exponent 
is the scalar product in M s . When Q has a probability density, the characteristic function is just its Fourier transform. 
The derivatives of C at the origin (if they happen to exist) are the moments of Q (if they happen to exist). Very 
helpful for our purpose is that it is easy to express the characteristic function for a scaled variable, say p,Q, with a 
fixed factor p,. Then we just have C m q(A) = Cq(//A). In our case Q — Q(t) will be the position after t time steps of 
a quantum walk W, starting from some initial state po, i.e., 

C Q{t) (\)=trp W t (e^). (9) 

When the walk is controlled by an external Markov process, we also have to take the expectation of the right hand 
side with respect to the stationary distribution of the control process. Now we want to look at a scaled position 
distribution. For example, in ballistic scaling we get 

C Q(t)/t (A)=trp W t ( e ^/ t ). (10) 

In the limit t — > oo the unitary operator exp(iA ■ Q/t) approaches the identity, which is invariant under W. On 
the other hand, we act on this nearly invariant element with a high power of W. The basic idea of our asymptotic 
evaluation is to look instead at the high powers of a slightly modified operator W e , defined by 

W E (X) = w(Xe leXQ y- i£XQ (11) 

where e = 1/t, or e = l/\/t for diffusive scaling, is now a small parameter. W £ is similar to W via the invertible 
linear (but quite non-positive) operator X H ► exp(ieA • Q). This means that 

w *( e feA-Q) _ w*(l)e teA ' Q (12) 

for all t. For initial states po supported on a finite region the exponential factor on the right hand side will be close 



to the identity in the scalings we consider, hence can be neglected when substituting this expression into ( 10 ) 



A crucial observation is that although exp(ieA • Q) is not a translation invariant operator, W e , like W commutes 



with translations. This is because if we apply a translation by x £ Z s to ( 11 ), we get two phase factors exp(±ieA ■ x 
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which cancel. In particular, if we apply W £ to a translation invariant operator like 1 we again get a translation 
invariant operator. However, we will consider W and W £ as maps on the space of translation invariant operators, 
i.e. as multiplication operators in momentum space, and apply perturbation theory in those subspaces. Now, the 
restriction of W £ to the translation invariant operators is not similar to the restriction of W. This is because the 
similarity transform on the whole space, i.e. right multiplication by exp(±ieA • Q) does not respect translation 
invariance. Hence, the eigenvalue of W £ which goes to 1 as s — > may differ from 1. Indeed, the perturbation theory 
of this eigenvalue is the core of our method. 

An operator A £ B(H) is translation invariant iff it is a function of momentum, i.e. (J c Aip)(p) = A(p)(Fip){p). 

Under the action of W or W £ such an A is transformed into a function A'(p). In general, A'(p) might depend on 
values A{p') at points p' ^ p. This is where Assumption n] comes into play: if there are no momentum transfers, 
i.e. each Kraus operator is itself a function of p, then (Wi)(p) = W(p)A(p) for a suitable operator W(p) on B{K) 
depending on p. Indeed, from 

(WA)(p) =Y,K a (pTA(p)K a (p) (13) 



and (11 1 we get 

(W E A)(p) = £ K a (p)*A(p)K a (p + Xe). (14) 

a 

The finite range condition makes each K a a trigonometric polynomial, so this operator is an analytic function of e and 
we can apply perturbation theory [3T]. Let us denote the Jordan decomposition of the above operator [HI Sect.I§5.4] 
by 

(W E A)(p) = ^(/^P^+D^))^)), (15) 

i 

where PiPj = (^jPj are the eigenprojections, and Dj eigennilpotent operators with D^P^ = PjD^ = <5j,-Dj for the 
eigenvalue (J>i(s). Let us assume for the moment that the eigenvalue Ho(0) = 1 of the unperturbed operator is simple, 
so that the unique eigenvector is 1, and the other eigenvalues satisfy |/Xj(0)| < 1 for i ^ 0. (This will be the standing 



assumption in Sect. VA but not in Sect. IV A I. Then Dq = 0, and we get 



Wl(l) = mo (e)*Po (£)(!) + •••, (16) 

where the dots stand for terms with i ^ 0. Since 1 is in the eigenspace for /i , and the e-dependent operators can 
be chosen analytic, all these contributions vanish as e — > , and Po(e)l — > I. Note that the crucial point here is the 



assumption |/i,(0)| < 1 for i ^ 0, cf. the discussion in Sect. VC So everything depends on the eigenvalue term (J,q(s) t . 



For ballistic scaling, i.e. e — 1/t and fJ,o{e) — 1 + iv ■ Xe + 0(e ), for some vector v g M s , we find 

Mo(e) ( =(l + ^+0(r 2 )) f -^e to -\ (17) 

Hence, the probability distribution of Q/t converges to a point measure at a deterministic (but possibly p-dependent) 
velocity v. A more detailed study of this case is given in Sect. |IV A"| and Sect. |V A[ where we also include an external 
control process. Consider on the other hand the special case v = 0, then the leading order contribution to {1q is of the 
form no(e) = 1 — \s 2 X ■ M ■ X + 0(e 3 ) for some matrix M. Then, in diffusive scaling s = t -1 / 2 we get 

This is the characteristic function of a Gaussian with covariance matrix M. Hence, the asymptotic dis tribu tion of 



Q/yi is Gaussian. A closer analysis of this case, again including control processes, will be given in Sect. VB 



IV. UNITARY QUANTUM WALKS 
A. Ballistic order 



The unitary case has been the subject of various papers [TJ 6; L -8, 10, ;42, 43]. Here we allow a general walk, as given 
by a unitary matrix W(p) (pi). As described in the previous section, we need to study the perturbation theory of the 
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eigenvalue 1 of the family of operators W e 



W e (A) = W( P )*AW{p + Xe), 



(19) 



where we treat p as a fixed parameter. Clearly, 1 is an eigenvector of eigenvalue 1 for this operator, but the eigenvalue 
1 is actually quite degenerate: any operator X commuting with W(p) is also in this eigenspace. The eigenspace is 
thus at least dim/C-dimensional, but if some of the eigenvalues of W(p) are degenerate at p, the degeneracy can be 
even higher. So let 



w{p) = Y J ^ k(p) Pk{p) 



(20) 



be the spectral resolution of W(p). In the case of degeneracies, i.e. when several 0Jk(p) coincide, this is not unique. 
But the perturbation theory of the one-parameter analytic family e n- W(p + eX) tells us that we can choose the 
operators Pk(p + eX) such that in the neighborhood of e = they and the corresponding eigenvalues u>(p + eX) depend 
analytically on e. In the sequel we assume such a choice has been made in (20). 

Not surprisingly, this leads to an analytic perturbation expression for W e . Indeed, let {P k }f = i an d {Rt\t=i be 
families of orthogonal projections in a Hilbert space K,. Then we can consider the operators E k g(X) = P k XRe on B(JC) 
one easily checks that each Em is hermitian with respect to the Hilbert Schmidt scalar product (Y\X) = tr(Y*X), 
and the E k t are themselves a family of orthogonal projections. Now setting Pf. = Pk{p) and Rk = Pk(p + eA) we find 



W £ = J2e lMp+eX) - u ' k< - p)) E k f, 



(21) 



kl 



which is clearly a spectral decomposition in terms of eigenvalues and eigenprojections, which are all analytic in e. 
Therefore, the expression 



w*(i) = J2 eitMp+£X) ^ k(p)) p k(p) p dp + 



(22) 



is correct to all orders. As e — > 0, we have Pi{p + eX) — > Pe(p), and since Pk{p) and Pi(p) are orthogonal, only the 
terms with k — £ survive in the limit. Moreover, with ballistic scaling e = l/ t the exponent converges to the derivative 
of uj(, which exists even at degeneracy points, because we have chosen (20) analytically. Hence 



hmW< /t (l) 



k 



exp^i 



dw k {p + eX) 
de 



(23) 



e=0 



Note, however, that the choice of projections P k at a degenerate point may well depend on the direction A, in which 
p is varied. Moreover, the derivatives of W(j>), compressed to the degenerate eigenspace need not commute, so there 
is no analytic choice of branches u>k- We call a point p a regular momentum for W(p) if in (20) we can choose P k 
and Wfc to be analytic functions of the vector near p. Of course, when the eigenvalues of W(p) are all non-degenerate, 
p is regular, and this will almost always be the case. 

For all regular p, we can write (23) as the exponential of the operator iX ■ V(p), where V is the p-dependent vector 
operator with components 



(24) 



This is the operator of group velocity. Note that, for all regular p, all of its components commute with W(p), since 
they are linear combinations of eigenprojections of W(p). Therefore, the components of V are jointly measurable in 
the sense of standard quantum mechanics. In any initial state p this gives a probability measure on velocity space. 
This measure is the asymptotic position distribution starting from p. 

Theorem 4 Let p \— > W(p) be the defining unitary of a quantum walk on M. s . Let Q(t) denote the position observable, 
evolved for t steps. Suppose that almost all p are regular, so the group velocity operator V(p) is defined almost 
everywhere. Then 



lim 

t— ^oo 



Q(t) 



v 



in the sense that for all bounded continuous functions f : R s — > C going to zero at infinity we have the weak operator 
limit of f(Q(t)/t, evaluated in the functional calculus, is f(V). This means, for any initial state p the distribution of 
Q(t)/t goes weakly (in the sense of probability measures) to the distribution of V in p. 
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Let us conclude this section with two short examples of unitary quantum walks without temporal disorder. We 
start with the generic one dimensional quantum walk with a two dimensional coin space, given by a shift operation S 
that moves the walker left or right depending on the internal degree of freedom and some SU(2) coin operation. In 
that case the walk operator in momentum space can be written in the familiar coin and shift decomposition 

W(n\ - <7 • <? - f ^(a) 6 '' 3 sin(a)e^ \ / \ , . 

with some unitary matrix C. By diagonalizing this matrix we get the p dependent phases u>± of the eigenvalues of 
the walk operator and the group velocities 



w±{p) = ± arccos(cos(p + a) cos(/3)) 
duj± , cos(/3) sin(p + a) 



dp 



± 



\J\ — cos(/3) 2 sm(p + a) 5 



In the case of the well known Hadamard walk, where C is chosen to be the Hadamard matrix, the dispersion 
relations are given by lu± (p) — | ± arccos (j^f^j , where the ^ compensates a factor of i due to our parametrization of 

the 2x2 unitaries. In Fig. 3(a) a plot of the p dependence of both eigenvalue bra nches is shown. From these dispersion 
relations we can compute the group velocities d p u>±(p), which are given in Fig. 3(b) Looking at the asymptotic t" 1 
scaled position distribution P(ir) of the initial state ip(p) = (q)j we see that at the points the density diverges 
(Fig. 3(c) Comparing this with the graph of the group velocities we find that this corresponds to points, where the 
group velocity has an extremal point. We call these points caustics and they are precisely the points causing the 
peaks in the asymptotic distributions. 




FIG. 3. (color online) For the quantum walk with Hadamard coin, the figures show (a) the dispersion relations ui±, (b) the 
velocities v±, and (c) the asymptotic probability distribution. The initial state is chosen to be i/>(p) = (J). The extremal points 
of the functions v±(p), hence the inflection points of w±, are responsible for the peaks in (c) at ±1/V2. 



Such caustics are also a generic behavior of quantum walks on Z 2 . In accordance with the one dimensional case we 
define a caustic a bit more formally as a point in momentum space where the Jacobi matrix of the group velocity or 
equivalently the Hessian of ui(pi,p2) is singular. This in turn implies that the density of the i -1 scaled asymptotic 
position distribution will diverge at such points, cf. Corollary [7J 

As an example we will study a quantum walk on Z 2 with a two dimensional coin space. The shifting in the p\ and 
Pi direction is done separately and in between two unitary coin operations are alternated. The overall walk operator 
is therefore given by 

W(p 1 ,p 2 ) = U 2 -S 2 -U 1 -S 1 =U 2 l e Q e-i^J-Oi-f o e -i Pl )- (26) 
The resulting dispersion relations are given by 



w ±(Pi,P2) = ± arccos(cos(pi + p 2 +6\ + 2 ) cos(</>i) cos(^ 2 ) 
- cos(pi -P2-X1+ X2) sin(0i) sin(0 2 ) , 



(27) 
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where greek lower case letters correspond the parametrization of a 2 x 2 unitary matrix as given in ( [25] ) and the 
subscripts refer to the unitaries Ui. From this relation we could also calculate the velocities and the points in the 
asymptotic distribution where caustics can be observed. We will do this for an explicit example, see Fig. [4] where 
the parameters of the two unitaries are chosen to be 6\ = —6 2 = §, Xi = — X2 = f , 4>i = f an d <fe = § 

The resulting band structure is depicted in Fig. 4(a)| In principle, the bands could cross or intersect each other 



at a line or a single point. In this example, however, we can observe the generic behavior of an avoided crossing of 
the two branches of the dispersion relation. The red lines in the graph indicate momentum values for which we find 
caustics, i.e. where the determinant of the Jacobi matrix of uj(pi 1 p2) vanishes. The second graph (Fig. |4(b)[ ) shows 






(a)Dispersion relation 



(b)Probability density of t)+ 



(c)Probability density of v+ 



FIG. 4. (color online) For the two-dimensional walk (261 the plots show (a) the dispersion relation ui±, (b) a contour plot, and 
(c) a 3D plot of the contribution of the cj+-branch to the asymptotic position density, for a particle starting at the origin. The 
red lines in (a) are curves of vanishing curvature. At these points the velocity density (i.e., the inverse of the Jacobian of the 
transformation p i-> v+(p)) diverges. This produces the enclosing red line in (b), and infinitely high values in (c). For more 
complicated walks such lines also appear in the interior of the velocity region. 



a contour plot of the possible pairs of velocities (d Pl u), d P2 u) in the upper branch of the dispersion relation uj + . The 
red line at the border of the possible velocities, where the value of the determinant drops to zero, corresponds to the 



red lines in the pi,p2~dependent graph in Fig. 4(a) of the dispersion relations. As in the one dimensional case, these 
are the points where the probability density of the asymptotic position in ballistic scaling will exhibit peaks. In the 
last graph (Fig. 4(c) I the value of | det(9pW + )|~ 1 is plotted dependent on the velocities. 



B. Higher orders 



We have evaluated asymptotic characteristic functions by employing a perturbation series. It is therefore natural 
to ask, whether one cannot use the higher terms in the series to get better approximations to the position density 
for large but not infinite t. Of course, this is possible, but one has to be careful in the interpretation of the results. 
The main problem is that the partial sums of an expansion of the characteristic function in powers of 1/t is not a 
characteristic function of any probability measure. Indeed, in the expansion we use we typically get combinations of 
\/t, so the higher orders of the expansion will be polynomials in A times an oscillating factor. This is clearly not 
integrable, so the inverse Fourier transform to get the probability density is ill-defined, and gives, at best, a rather 
singular distribution. However, if we only look for the expectations of sufficiently smooth functions of velocity, say 
f(Q(t)/t), we get the integral of the characteristic function with the Fourier transform of / which decays rapidly 
enough to absorb all polynomial factors. Thus for a fixed smooth test function the expansion makes sense. Another 
way to look at this is to multiply the expanded characteristic function with a suitable cutoff-function (enforcing 
sufficient decay in A) before transforming back to velocity space, resulting in a smoothed out probability density. 
Then "removing the cutoff" and the series expansion do not commute, and the choice of cutoff is effectively the choice 
of a smooth family of test functions. 

We can go back directly to Eq. (22 ), which is correct to all orders. That is, the characteristic function of the position 
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distribution in ballistic scaling at time t is 

t 

1/ 



C t (X) =trpWi /t (I)e iA -«/' 



= /"dp ^ e ^(^(P + A A)-^b)) tr ( p (p + AA,p)P fe (p)P £ (p + AA)) (28) 

Here, as in the previous subsection p(pi,P2) denotes the integral kernel of the initial density, and the trace is over the 
internal degrees of freedom. In leading order we could neglect the shift by X/t in Pi, so only terms with k = £ remain. 
Looking now at the first order term resulting from the expansion of Pi(p + X/t) and k ^ £ we find an oscillatory 
integral with a regular integrand and rapidly oscillating exponential exp it(w(,{jp) — iuj k {p)). Assuming that this phase 
is not constant on sets of positive measure (as a function of p), we conclude that the integral goes to zero, so that 
with the factor 1/t from the expansion of Pi such terms are o(i -1 ) and can be neglected. Of course, in higher order 
corrections one will have to extract the leading orders of the oscillatory integral by a stationary phase analysis. 
For the expansion to first order we need the expansion of the dispersion relation to second order: 

u k (.P + X/t) = uj k (p) + - ■ v k (p) + 7^u"(p, X) + o(t- 2 ), (29) 

where uj" is a quadratic from in A containing the Hessian of the branch uj k of the dispersion relation. This approx- 
imation eliminates one of the most prominent features of the finite t probability distributions, namely their rapid 
oscillations. Indeed, from inspecting such distributions it is clear that these oscillations have a frequency of order 1/t 
in ballistic scaling, i.e., they are really at the scale of the underlying lattice. This is reflected in the exact formula 



(281 by the fact that all expressions are 27r-periodic in X/t. So the Fourier transform of Ct is a sum of (^-functions 



at the lattice points. The approximation (29) destroys this feature, resulting in a rather smooth function inside the 
allowed region of velocities. The Hessian can be determined from standard second order perturbation theory of W(p), 
but we found it convenient to eliminate it by partial integration of the 1/i-term containing uj". Since the integrand is 
periodic in p, this gives no boundary terms, and the resulting differentiations of the trace can be combined with the 
other Taylor expansions to give the first order term 

C t (A) = C oc (A) + ^ J dp^e^Phr( K (Q-X)p + p(Q-X))(p,p)+ t [P k (p),X-P k ( P )})+o(t- 1 )^ (30) 

In the example of a Hadamard walk one can substitute the variable v(p) and find a Bessel function for the first order 
term. However, due to the factor A this is not integrable, and the inverse Fourier transform has to be taken in the 
distributional sense, more specifically under the integration with smooth test functions whose support in velocity 
space stays away from the caustics. This results in the expression 

Pt( U ) = -7j ^ + 7 ■ = 3 + fOT M < ( 31 ) 

tt(1 - u)Vl - 2u 2 t 7r ^i_2u 2 V2 
for the Hadamard walk starting at the origin in state ipo = (1, 0). This is shown for t = 10 in Fig. [5j 

V. TIME DECOHERENT WALKS 

A. Ballistic determinism for decoherent processes 

We now extend the analysis to decoherent models, as described in Sect. |II A| We will continue to work under 
Assumption [l] i.e. excluding momentum transfer. We can thus still apply the theory outlined in Sect. |III[ and find 



that the asymptotic distribution of position is given, as in ( 14 1 by an operator depending on momentum. This is 
also true if we include a control space with a Markov chain dynamics. The observable space we work on is hence the 
tensor product of C°°(T), i.e. the control space variables with the quantum observables B(fC). We will fix the value 
of p, and everything will depend on this parameter, but this dependence will, for the moment be the only thing left 
of the translation degree of freedom. 

As before, we will consider the observables on A <G C°° ®B{JC) as measurable functions 7 1— > A(j) 6 B(JC). The two 
ingredients of each time step will then be on the one hand the 7-dependent quantum operation 

(V4)( 7 ) = V 7 (A( 7 )) = £x; Q (p)A( 7 )ir 7 , a (p), (32) 
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P(x) 




FIG. 5. (color online) Correction of order 1/t to the asymptotic position distribution for the Hadamard walk from the initial 
state Vo = (liO) located at the origin. The polygon connects the exact values for n — 10. The asymptotic distribution 
(red/dashed) overestimates the left peak and underestimates the right peak. The correction after ( |31[ ) for the same t is shown 
by the green/solid curve. 



where the Kraus operators are normalized such that the sum over a with A( 7 ) = 1 gives 1. The second step is the 
update of the control parameters by the Markovian evolution M, i.e. 



(MA) ( 7 ) = J m 7 (d7 ? )A( ?7 ). (33) 
The full evolution is then given by W = VM or, written out more explicitly: 

(WA)( 7 ) = J m 7 (dr 7 )V 7 (A(r ? )), (34) 



with the momentum dependence implicit in the dependence of V 7 according to ( 32 ) . 

Obviously, this form contains the unitary case, e.g., when V 7 is the same walk unitary for all 7, and the control 
process is irrelevant for the walking particle. We are now interested in the opposite end, i.e. the case of generic 
randomness. Therefore we will assume the following 



Assumption 2 For almost all p, the eigenvalue p,Q = 1 o/W as an operator on C°°(T) ®B{fC) is simple and isolated 
and for i =^ the eigenvalues satisfy \pi\ < 1. Moreover, there exists an invariant faithful state for W, i.e. a density 
operator p with nonzero eigenvalues on K. such that tr p V 7 (A) = tr p A for all 7 and all A. 

An example of commuting Kraus operators, where the non-degeneracy condition of the eigenvalue 1 is violated 
given in |VI E| To see that this assumption is not overly restrictive we consider the following scenario. 

Proposition 5 Assumption^ is valid if the quantum walk W has the following properties 

1. r is finite, so the transition probabilities can be written as a matrix m 7 (ry). For some power of this matrix all 
entries are strictly positive. 

2. There is a density operator p with nonzero eigenvalues on K, such that tr p V 7 (A) = tr p A for all 7 and all A. 

3. For almost all p, and some r € N the set of operators K liai (p) ■ ■ ■ K lr a r (p) is irreducible, i.e. only multiples of 
1 commute with all of them, and its linear span contains the identity. 

The last condition seems tricky to check, but it is generically satisfied. In fact, the operator products usually span 
the whole space of matrices for relatively small r. 

Proof Since £°°(r) Cg)i3(/C) is now finite dimensional, any simple eigenvalue is isolated, so we only have to prove that 1 
is a simple eigenvalue, and the only one on the unit circle. The transitivity assumption on the transition probabilities 
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guarantees that there is a unique, strictly positive, invariant probability distribution m on T. Hence we have an 
invariant state for W, namely m (g> p, written out as 

(m®p)(A) = J2 m y A A (i))- 

7 

We first show the simplicity of 1, and since linearly independent eigenvectors of W will be independent eigenvectors 
of W" we may do this for some power of W. We choose some multiple of r, say rn, chosen sufficiently large so 
that all entries of the nr-step transition matrix are positive. Then condition 3 in the Proposition is also satisfied for 
rn, because the identity lies in the span of the Kraus operators for W. Hence, for this step we can simplify the 
assumptions to m^{rj) > for all 7,77 and the irreducibility of {K ja }. 

The basic technique for the proof is the decomposition 

W(A*A) - W(A)*W(A) = v(m{A*A) - M(A)*M(A)) + V(B*B) - V(B)*V(B), (35) 

where B = M(j4). Both terms are positive by the "2-positivity inequality" [3S] for channels. But when we evaluate 
for W(A) = A in an invariant state of W the left hand side becomes zero. This will provide a lot of information from 
the vanishing of sums of positive terms on the right. Explicitly we get 

(M(A*A) - M(A)*M(A)) (7) = \ m 7 (77)m 7 ( x ) (A(r,) - A( x )) * (A{ V ) - A( X j) , 

riiX 

where we used that m 7 (r;) = 1 for all 7. Applying the invariant state m (g> p , and using the invariance condition 
for each V 7 , we find 

m7 m 7 (v)**h (x)trp(A(ri)-A(x)y{A{r])-A(x)) =0. 

Since the probabilities m and m are all strictly positive, we find that each summand vanishes. Because p has no zero 
eigenvalue, this also implies that (A(i]) — A(x)) (A(n) — A(x)) = 0, and hence A(rf) = A{\) for all T],X- Hence we 
can set A{^f) = A for all 7. This makes the first term in ( [35] ) vanish for every A, and leads to B = M.(A) = A in the 
second term. Moreover, from p4| ) we see that W(A) = A just means that V 7 (^4) = A for all 7. 
Now consider, for each one of the operators V 7 , which has Kraus operators K~ a , the expression 

J2lAK ia }*[A,K ia ] =Y 1 {A*A)-\ 1 (A)*A-A*\ 1 {A) + A*A (36) 

a 

If A is invariant, this reduces to V 7 (A*A) — A* A, which is clearly zero under the common invariant state p. Hence 
the expectation of the positive terms on the left hand side must vanish also, and since p has no zero eigenvalues each 
[A, Kj a ] = for all a and 7. By assumption this implies that A is a multiple of the identity. 

It remains to be shown that there are no further eigenvalues on the unit circle. Suppose to the contrary that 
W(X) = w!/0 for some u ^ 1 with \uj\ = 1. Then the operator W(X*X) - W(X)*W(X) = W(X*X) - X*X, 
which is positive by the 2-positivity, has vanishing expectation in the faithful invariant state, which implies that it is 
zero. But then X*X is a fixed point of W, and we have already seen that this implies that X*X is a multiple of the 
identity. Since it cannot be zero, we can normalize X so that it becomes unitary. But 2-positivity also implies that if 
W(X*X) - W(X)*W(X) = we must also have W(Y*X) - W(Y)*W(X) = 0, for all Y. Otherwise, the inequality 
could not be valid for linear combinations of Y and X. But then, by induction on n, we find that W(X n ) = ui n X n . 
In other words, all powers of ui are eigenvalues. Since the dimension of the space is finite, this means that w must be 
a root of unity, say uf 1 — 1. But then X is also an eigenvector of W™ with eigenvalue 1. Since our arguments for the 
simplicity of 1 also apply to powers of W, this implies X = 1, and a contradiction to w ^ 1. ■ 

Now we can apply the ideas of Sect. |IV A| to extract the ballistic spreading of the position distribution. The only 
difference is that we can now use non-degenerate perturbation theory, so everything is much simpler. Consider again 
some parameter A as the argument in the characteristic function, and the operator W e , acting on A E C°° (T) ® B(H) . 
We define operators V xe via 

(W E A)(j) = [ m y (d v ) Y.KyaipTAi^K^ip + eX) = f m y (d V ) V 7 , e (A(r?)) . (37) 
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Let W £ A e = p e A e denote the branch of eigenvectors and eigenvalues with po — 1. Since the eigenvector is only 
determined up to a factor we are free to choose a normalization so that 

m®p(A e ) = l. (38) 

Expanding these objects to second order, we get 

v 7 , e = v 7 + e v; + ^v;' + o( £ 3 ) 

e 2 

p £ = 1 + ep! + — + 0(e 3 ) (39) 
A e = 1 + eA' + —A" + 0(e 3 ) . 
For now we will only use the first order, in which the eigenvalue equation reads 

W(A')( 7 ) - A'( 7 ) = ,/I - J m 7 (dry) v;(I) = //I - V' y (t) 



ti'l-J2 K -y» (p) * ^K ja (p + e\) 



(40) 



e=0 



Where we evaluated the sum over rj, on which the integrand does not depend. Taking the expectation with respect 
to the invariant state m <g> p makes the left hand side zero, and leaves an explicit equation for p! , namely 

p = i\-v(p) with the vector 

v(p) =-i m(d 7 )^trpA- 7Q (p)*Vif 7Q (p). (41) 

J a 

This expression is real, because for all 7 and all X € B(JC), the function 

y j tr p K ja (p)*XK ia (p) = trpX 

a 

is constant, and hence has zero gradient. Summarizing we get the following result. 

Proposition 6 Given Assumptio n^]] and Assumption the asymptotic distribution of Q(t)/t is determined by the 
distribution of v(p) as defined in Mjy - 

Proof According to the preceding reasoning the characteristic function of the random variable Q(t)/t converges for 
t — > 00 to f dpe lX ""( p ) ti p(p) and the asymptotic distribution of Q(t)/t is given by inverse Fourier transform of this 
function. ■ 

Two features of this formula are remarkable. Firstly, v(p) is a scalar, i.e. a multiple of the identity with respect to the 
internal degrees of freedom from K,. Therefore, in contrast to the unitary case discussed in Sect. |IVX} the asymptotic 
distribution is independent of the initial state of the coin. Secondly, the transition probabilities of the control process 
only enter through the invariant probability distribution m. This means that we get the same result for a Bernoulli 
process, in which we take independent coins in successive steps with probability distribution m. The control space 



therefore just contributes another index to the Kraus decomposition of a one-step channel, cf. Sect. VD We will see, 
however, that the transition probabilities are not irrelevant for next order, see Sect. |VI A| 

The following corollary gives a closed formula for the asymptotic distribution of Q(t)/t if the function v(p) is such 
that there is a decomposition of momentum space into a finite number of open sets on which the function v(p) is 
invertible and a boundary which coincides with the set of caustic points of the quantum walk. 

Corollary 7 Let the characteristic function of the asymptotic distribution ofQ(t)/t be given by 

lim C Q(t)/t (X) = [dp e iX <rtp(p) 

( — 7T,7r] s 



17 



and suppose there exist finitely many disjoint open sets Rj, j G J on which v(p) is invertible and R = ORj is the set 
of all points such that the determinant of the Jacobi matrix dv(p)/dp is non-zero. Then, the asymptotic probability 
distribution ofQ(t)/t is given by 

-l 

{x) P {vf\x)) ,xeR s , (42) 
with Vj the restriction of v to Rj and the set of caustic point coincides with the complement of R in (— 7r,7r] s . 



P(x) 



E 



det U 



Proof The asymptotic position distribution P(x) can be obtained from the characteristic function C(A) by applying 
the inverse Fourier transform 



P(x) = Jdxjdp e a «^p{p) . 

W s ( — 7T,7r] s 

Using the integral representation of the Dirac distribution S(x) = t^jf Jgs dX e lX ' x we get 



P ( x ) = E dp S ( V (P) ~ x )p(p) ■ 



Now, formula (42 1 follows from the fact that v(p) is invertible on the Rj and a substitution of the integration variable p 
by v. The complement of R in (— 7r, 7r] s is given by the set of points where the Jacobi matrix dv(p)/dp, i.e. the Hessian 
of the dispersion relation, is singular. This is the set of caustic points, and the asymptotic probability distribution 
diverges at the corresponding points. I 



B. Markov controlled walks in diffusive scaling 

We now go on to study the asymptotic position distribution on the \ft scale. If Q{t)/t — > v goes to a sharp value 
the quantity of interest will be (Q(t) — vt)/y/i, and we can hope that this has a well-defined limiting distribution. 
Indeed for this to make sense we need that the velocity is sharply defined. Otherwise, we would just see the limiting 
distribution of v scaled to larger and larger variance. Therefore, we need to restrict to a setting, in which v{p) is 
automatically constant. So from now on we assume that we have a Markov controlled random coin as described in 
the introduction: 

Assumption 3 For every 7, the operation V 7 is unitarily implemented. That is, for each 7 there is only one Kraus 
operator K y — Kyi, which is a unitary walk operator, so that V 7 (A) = K*AK~f. 

An example violating this assumption is given in section |VI C| The following proposition says that this indeed 
implies the required constancy of the ballistic velocity. 

Proposition 8 Assumption^and Assumption^together imply that the function v(p) , defined in is independent 
of p. Moreover, 

' 1 and (43) 



dim/C 



■-/■"wSa?- (44) 



Proof Obviously, the density operator p oc 1 is an invariant state for each V 7 . By definition of the index, det (p) 
cexp(zp • ind K^). Hence 



iX ■ ind -fC-, 



det Ky(p) de 



det Ky (p + eX) 



e=0 



= — det(iY 7 (p)*if 7 (p + eA)) 



e=0 



tr (Ky (p)* — Ky (p + eA)) = dim K tr p V 7 (1) , 
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and the formula follows by dividing this equation by dim/C and summing with respect to m. ■ 

Note that for unitary quantum walks W the index ind W is given by the trace of the group velocity operator V, 
see Sect. |IV A Therefore, if we naively apply (44 1 to the unitary case, although Assumption [2] is violated, we find 
v(p) = tr V(p). Of course, the asymptotic distribution of a unitary quantum walk W is not determined by v(p) but 
by the whole operator V(p). 

When the ballistic order is completely defined by a deterministic velocity, i.e. Q(t)/t —> vl, then we can subtract 
the ballistic motion and look at how the probability distribution develops around it. In other words, we look at the 
deviation operator 



D(t) = ^=(Q(t)-vt). 
We will compute the limit of the characteristic function of D, i.e. 



lim Cd( ( )(A) = lim e 



— i\-v yi 



a 



lim e" 

t—too 



iX-v \/t 



Q(t) 



dp 



A 

Vt 



I(d 7 )tr (ptpJW^I)^) 



dptv(p(p)) lim e~ iX v/Vi (l 



.X-v n" 



0(t" 3 / 2 ) 



dptr(p(p))exp ( -((j," + (\. v ) 2 ) 



(45) 



(46) 



Here, at the first equality we substituted ( |45|, a t the second introduced W from ( |11[ ) in the form (37), and at the 
third introduced the perturbation expansion (|39|) with e = \j\fi. At the last equality we used again the asymptotic 
formula lim t _ i . 00 (l + x/t + o(i -1 ))* = exp(:c), which follows immediately from the Taylor expansion of the logarithm. 
In the application above there is a cancelation of large phases, so the formula must be applied with care, in the form 



lim e 

t— too 



ia\ft 



= lim I 1 

t—>oo 



Vt 

b 
t 



1 + ^ + - + 0(t- 3 / 2 ) = lim e -'VVt(i + ^ + - + 0(t- 3 / 2 )) 



2t 



+ o(t- L ) =exp [b 



a 
~2 



(47) 



Note that, in contrast to p! — iv ■ A, the second order perturbation coefficient /z" will depend on p. Moreover, since 
the perturbation is proportional to eA, the second order perturbation coefficient fjf' must be homogeneous quadratic 
in A. Hence, according to (46 ), the limiting distribution of D(t) is a mixture of Gaussians with p-dependent covariance 
matrix s(p) such that 



A • s(p) ■ A = y% ifc (p)AjAfc = -fi"(p) - (A-u) 2 
jk 



(48) 



For the moment, we fix A, and concentrate on computing the expression (48) for any set of given data and verifying 
its positivity. 



Evaluating the eigenvalue equation ~W £ A e = fi e A s to first and second order in e, we now get 
W(A')(7) - A' (7) = iX-vl - v;(I) 

W(^0(7)-^ / (7)=M // 1 + 2m'A / ( 7 )- / m 7 (dr ? )(v;'(l)+2V;(^(7 7 ))) 



(49) 
(50) 



The first line is just a repetition of (40). From the second line we extract \J' by taking the expectation with respect 
to the invariant state, using also that by the convention (38) the second term on the right has zero expectation: 



M "= / m(d 7 ) tr(pV;'(l))+2 mid^m^dr,) tx(p V^A'fa))) 



(51) 



In this equation, the first order perturbation A' of the eigenvector must be extracted from ( 49 1 . Indeed, this is 
uniquely possible: By Assumption [2] the simple eigenvalue 1 of W is isolated, so the rank of W — id is exactly one less 
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than maximal. The kernel is explicitly known: W — id annihilates precisely the multiples of the identity and maps 
onto the elements with vanishing expectation under migip. By definition of v, the right hand side hence lies in the 
range of W — id, so there is a unique solution A' satisfying the normalization condition ((38]). Note that since the 



solution is unique, and is skew hermitian, so is A'. Moreover, since the right hand side of (49) is linear in A, so is 
A'. We stress that Equations (49)-(51 ) are also valid if Assumption [3] is violated and for ballistic scaling e = 1/t with 
non-constant v(p) they can also be used for a finer analysis of the large time behavior, incorporating second order 
perturbation theory. From now on we will adopt Assumption g so that V;(l) = K*K' T Since each K~ i is unitary, 
we can use this to express the derivatives of K 1 by A': 



K' y = K 7 (A'( 7 ) 



W(A'){j) + iX-vl) 



(52) 



We come back to the determination of p" from |51J). We can eliminate the second derivative in the first term by 
differentiating (41 ) with respect to p, or more precisely, by setting p = p + eX and differentiating with respect to e. 
This gives 



/d 
m(d 7 ) — tx(pK*(j> + eX)K' 7 (p + eA)) 

m(d7) tv{pK*'K') + -i /m(d 7 ) ti(p\''(Tj) 



and hence, using (52) and the unitarity of K 1 



\{di) tr(pV;'(l)) 



m(dj) tr(p \A'(j) - W{A!)(i) + iX-vlf), 



(53) 



(54) 



where, for an operator X, we use the abbreviation \ X\ 2 = X*X. The second term in (51) can also be simplified by 
eliminating the derivative in V^(X) = K*XK' 7 via (52). We get 



J m 7 (dr?) k;A'( v )K^ = W(A')(7)(A'(7) - W(A')(j) 



iX-vl) 



Note that in (51) we need the expectation of this expression in the invariant state m (g) p, just as (54) is such an 



expectation. Bringing together the various terms of (51), and using the skew hermiticity of A', W(A)', and iX-v, we 
find for p8|: 



p" - (X-v) 2 = = m ® p (\A' - W(A') + iX-vl\ 2 - 2W{A'){A' - W(A') + iX-vlj) - {X-vf 
= m(g>p (w{A'f - (A') 2 + [A', W(A')] - 2iX-vA'^ 



m(g>p(\A'\ 2 - \W(A')\ 2 



(55) 



Since m ® p is invariant under W , we can also write the first term as the expectation of ), so that this 

expression is non-negative by virtue of the Cauchy-Schwarz inequality for channels [35]. With A' a linear function of 
A, the above expression becomes a quadratic form in A, as claimed in (48). 



Intuitively, one would expect that the diffusion constant becomes very large if the considered quantum walk differs 
only little from a coherent quantum walk with ballistic behavior. For example, if the coin operations of a decoherent 
quantum walk would be all very similar or if the Markov process is such that it prefers one of the coins, we would 
expect reminiscences of a coherent quantum walk even for large times. This effect can be seen in the examples in 
Sect. VIA and VI B where we derive an explicit formula for the variance s(p), which diverges in the coherent limit 



of the the considered quantum walks. Similarly, if the Markov process converges to a deterministic Markov chain, 
i.e. the transition matrix approaches a permutation matrix, the diffusive order will diverge. This is because we 
could consider a full cycle of the permutation as one step of another quantum walk, which is also unitary, leading to 
ballistic transport. If on the other hand the transition rates of the Markov process are very small, i.e. there are long 
subsequences of coherent evolution, the overall evolution up to time step t would be very similar to the average of 
a number of coherent quantum walks each evolved for t time steps. If some of these coherent quantum walks show 
ballistic behavior we would again expect that the diffusion constant diverges. 
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C. Higher orders without ballistic determinism 



The general remarks about higher order expansions in Sect. IV B| apply also in this case. In this section we focus on 
just the first (i.e., 0(l/t)) correction to the ballistic scaling. When the ballistic velocity is independent of p, this is, in 
fact, best expressed by the diffusive scaling. However, almost all of the second order perturbation theory developed in 
the previous sections is independent of that Assumption [3j Therefore, we have already done most of the work needed 
to get the first correction to ballistic scaling in the general case. 

We adopt Assumption [T] and Assumption [2j but not nece ssarily Assumption [3] Let us neglect for a moment the 
Markov control, and consider the Jordan decomposition (151 for small e. The t th power of this operator is built from 
the terms 



(A*i(e)Pi(e)+D<(e) 



r=Q 



Here is the order of nilpotency of D-j(e), which is bounded by the algebraic multiplicity of Pi(s). When we choose 
the label i = for the key eigenvalue (i.e., /xo(0) = 1) we have < 1, for all i by Assumption [2j Hence 
expressions such as t n p l ~ r go to zero faster than any power of t~ n (n € N). Consequently, in any expansion in such 
powers, all terms but the one with i = can be neglected. Moreover, p is simple, so r = 0. Therefore, only the 
term P (£)^o(£)' remains. Since we were only interested in the leading order so far, it was enough to replace Po(e) 
by Po(0). However, for a systematic expansion we also have to expand the rank one projection Po(e) in powers of 
e. This will give Pq(e)(X) = A(e)tt(pB(e)X), where A(e) is the eigenvector of W e (p) and B{e) the corresponding 
eigenvector of the adjoint. Of course, eigenvectors are only defined up to a (possibly e-dependent) factor. We have 
already used the convention that trpA(e) = 1. This fixes the factor also for B(s), since we must have Po(e) 2 = Po, 
and hence tvpA(e)B(e) = 1. Hence, expanding B as in p9| we get trpA' = tipB' = 0, so that for an expansion to 
order t~ l we do not need B'. These ideas remain valid with Markov control, with the invariant state p taking 
the role of p . Then 

P (i)( 7 ) = l+^'( 7 )+o(<- 1 ), (56) 



t 



where A' is determined exactly as above. For the eigenvalue /xo(e)' = (1 



2t 2 



we first expand the logarithm, 



giving cxp(// + f (p" - (p') 2 ) + o(t r ). 



p (e) f = e 



With the p-dependent covariance matrix s (see (48)) we can thus write 
v(p) 



iX 



1 \ 

Tt A 



s(p) ■ X 



JX ■ v(p) i 



— A • s(p) 
2t yF ' 



or 1 ) 

a) +o(r 1 ) 



(57) 



From the point of view of a series expansion in inverse powers of t these two forms are equivalent. However, the first 
form is preferable, because it is the characteristic function of a Gaussian distribution and hence has a probabilistic 
interpretation. This has to be taken with a grain of salt, however: The proof of positivity of s{p) given in the previous 
section does depend on Assumption [3j Indeed, we will see in an example (Sect. VIC) that s(p) may be complex, and 
only represents a probability after integration of the whole expression with respect to p. 

Finally, we have to expand the second factor in (12). This correction contains some information about the initial 
position distribution, although contracted by a factor 1/t. When p(p\,p2) denotes the integral kernel of the initial 
density, the effect of the factor exp(iA • Q/t) under the trace is to shift the first argument of p. We therefore introduce 
the function 



C (A,p) = trp(p + X,p), 



(58) 



where the trace is over the internal degrees of freedom. This notation is to suggest that this is some kind of char- 
acteristic function of the initial distribution. This is not literally true, since we keep the momentum variable, so at 
best it is a phase space distribution function. Indeed, a slightly more symmetric version ti p(p + \/2,p— A/2) is the 
Fourier transform of the "position distribution at fixed momentum" according to the Wigner distribution function. 
Of course, this distribution is contracted by a factor 1/t in ballistic scaling. Together, we get 



C t (A) = dp e 



JX ■ v(p) 



l \ 

2i A 



s(p) - A/ 



(c (X/t,p) + \tvp{p,p) { m{di) A'( 7 ) 



o(r 



(59) 



Hcuristically, the interpretation of the first term is the distribution of the sum of two independent quantities for each 
p: a Gaussian centered at v(p), with decreasing variance, plus a scaled down version of the initial position distribution. 
This is then averaged over momentum. The second term does not allow such a simple interpretation. Two numerical 
examples are shown for quantum walks in one lattice dimension in Sect. VI A and VIC below. 
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D. Memoryless decoherence 



In many applications there is no memory in the control process: the walks V 7 are chosen independently in each 
step. The Markov process is then a Bernoulli process, and the probability m 7 (?j) to end up in 77 is the same from any 
state 7. Obviously, this probability is then also the invariant distribution, i.e., we have 

m 7 (7y) = 111(77) (60) 

for all 7,77. We can look at the resulting process in two ways: On the one hand, we could just specialize using ( pj| . 
which is the approach we will take below. On the other hand we could completely discard the control process, so 
there is only one value 7 = 0, say, and the corresponding V operation is V = fm((2 7 )V 7 . This would violate 
Assumption 3 but, of course, this case can nevertheless be fully analyzed. In fact, it simplifies the computation of the 
operators A' and the diffusive order fj," considerably. We can reduce the dimension of the system of linear equations 
that determines by a factor n, where n is the number of Kraus operators from which the control process can choose 
from. 



Substituting (60) into equation (51| for fx we get 

fi" = / m(d 7 ) tr(pV 7 (l)) + 2 / m(d 7 ) k(p%{ \ m(d V )A' (r,))) 



(61) 



So in order to determine //' we do not need to know the individual A' (77), but it suffices to know the value of the 
average 



X = A' 



m(drj)A' (ij). 



(62) 



Our aim is to set up an equation directly for this unknown matrix X. From (34) we find (WA')( 7 ) = V 7 (X), which 
turns (49) into a definition of A' (7) in terms of X 



A'(rr) = v 7 (x)-iA-«i + v;(i) 

The equation for X now follows by averaging: 

m(d 7 )V 7 (X) - X = iX-vl - I m(d 7 )V' (1) 



Eliminating the control process from (64) and (61 ) via the definition of Vq and X we get the equations 



and 



V (X) - X = i\-vl-V' (l) 
Ai" = tr(?V£(I))+2 tr(pV' (X)) 



(63) 

(64) 

(65) 
(66) 



VI. EXAMPLES 



The Hadamard walk with reflections 



The aim of this section is to derive the asymptotic distribution of Q(t)/y/i of a quantum walk where the control 
process is of Bernoulli type and in each time step it picks one of two unitary quantum walks. We will highlight the 
simplification of the formulas (49) and (50), accounting for general Markov processes, for the case of Bernoulli type 



decoherence by also calculating the diffusive order for general Markov processes, which then doubles the dimension 
of the computation. The unitary quantum walks, the control process chooses from, admit a decomposition into shift 
and coin operation. The shift operation is assumed to be the same for both walks, hence, the control process chooses 
only the coin operation for each time step. With probability 1 — e it will pick the Hadamard coin for each time step, 
which, on its own, would lead to ballistic behavior. In the remaining case the Pauli matrix o\ will be applied, this 
unitary quantum walk hinders the particle from moving at all, the particle will be reflected at each site. The Kraus 
operators of the quantum walk read 



Ki(p) = 



1-e 



- ip 



K 2 (p) = \/i 







- 1 p 



(67) 
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First, we note that the velocity v(p) is zero, which follows from Eq. (41) and V'(l) = K{K[ + K^K 2 — iXa 3 . 
Alternatively, one could also rephrase this Bernoulli type walk as a Markov controlled walk and then use formula ( 44 ) 
to get v(p) — 0. This means, in ballistic scaling the position distribution converges to a point measure at the origin. 
In order to determine the diffusive scaling we make directly use of Eq. ( 66 ) . It is easily seen that 

V"(l) = K\K'{ + KIK'i = -A 2 1 
and hence by (JsTJ) and K[K* = (1 — e)a^ together with K 2 K 2 = £0-3 

s(p)A 2 = -//' = A 2 - tr ( K * A ' K a) = A 2 - ((1 - e) tr (Ha 3 H*A') + e tr (aia^ A')) 



with the Pauli matrix a\ and the usual Hadamard matrix H. Using the rules H* a 3 H = o\ and <7\<j 3 o\ = — 03 we get 

s(p) = A 2 - i2\ (r(l - e)oi - ea 3 ) 
where we abbreviated at — 2 _1 tro^-A', which implies A' — o-i^i- Now, A' is determined by (65 1, which reads 

3 

^ Oi (K\oiK x + K 2 a l K 2 - tJi) = -iAcr 3 . 

i=0 

By calculating the Hilbert-Schmidt scalar product of this equation with the Pauli matrices c.; we get the following 
system of equations 

e cos(2p) - 1 sin(2p) (1 - e) cos(2p) A 
£sin(2p) — cos(2p) — 1 (1 — e) sin(2p) ■ I a 2 
l-e -l-e 

where we have chosen ao = in order to satisfy m ® p (A 1 ) = 0. The solution to these equations is a\ = a 3 = iX/2e 
and a2 — iAtan(p)/2e. This yields the following expression for the variance 











)■ 




- 













sb{p) 



l-e 

e 



By inverse Fourier transform of the characteristic function lim Cof^fA) = e 
position distribution in diffusive scaling 



-\ 2 (l-e)/(2e) 



we obtain the asymptotic 



P(aO 



2tt(1 - e) 



The 1 jt correction to the asymptotic position distribution in ballistic scaling can be inferred from ( 59 1 , we omit the 
computation and just give the resulting probability distribution (cf. Fig. [6]), which reads 



Pt{x,e) 



2-Kte 



1-eV 2(1 -e) 



exp 



tee 
2ll~e) 



(68) 



For a general Markov process M, which chooses from the two unitary quantum walks, the time evolution reads 

W(A)( 7 ) = Y, nhiriW^Ain)^ , 7 = 1,2 
77=1,2 

now with Kraus operators and transition matrix 

The Bernoulli control process corresponds to the choice mi = 1 — m 2 = 1 — e. Since K\(p) and K 2 (j>) are irreducible 
for almost all p Proposition [8] applies. Observing that mdK y = we get v(p) = and hence Q(t)/t converges to a 
point measure at zero. Before calculating the diffusive order of this quantum walk we will prove a useful lemma. 
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FIG. 6. (color online) For the quantum walk according to (67 1 the plot shows the 1/t correction to the asymptotic position 
distribution in ballistic scaling after t = 10 depending on the parameter e. 



Lemma 9 Let K 7 ,7 6 T be a finite collection of unitary walk operators on Z with two dimensional coin. Assume 
the K 7 admit a shift coin decomposition K-y = U 7 ■ S with 



S(p) 



e'P 
e~ ip 



and p independent unitaries U 1 such that Assumption^ is satisfied. Let the overall time evolution W be given by 
Markov process M that chooses from the K J7 i.e. 

W(A)( 7 ) = J2^(v)K;A(r 1 )K y . 
r/er 

Then v(p) = and the diffusion constant is determined by 

A 2 • s(p) = -A 2 - 2Ai y~] m 7 a 3 (7) , 

■yer 



where 03(7) = 2 1 tr(cr 3 A'(7) with A' according to {49) and the Pauli matrix 0-3 



Proof Since ind = for all 7 it follows from (44) that v(p) = 



According to (49) and (55), we need to solve the equation 

W(A')(7)-^'(7) = -V;(l) 

in order to determine the diffusive scaling of W. To begin with, we decompose the components A' (7) into Pauli 
matrices c,- , i = 0, 1, 2, 3 with ctq = I2 



i=0 

Next, it is easily seen that for arbitrary two dimensional unitaries U~ l we have K*K' 7 = iXa^ and hence 

2 

W(A')( 7 ) = (03(7) ~ *A)a 3 + ]T 04(7)0* . 



i=0 



Now, the variance s(p) can be computed from (55) 

A 2 • s( P ) = m®p | 2 - |lf(i')| 2 ) ^E^ tr (l^'WI 2 - WW)! 1 



7er 
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and since tr \A'(j)\ 2 = 2 • J^Uo Mt)| 2 [t follows 



A 2 • a(p) = m7(l«3(7)| 2 - I« 3 (t) - *A| 2 ) = -A 2 - 2Az £ m 7 a 3 (7) , 
which is real since -A' (7) is skew-hermitian. 

By a straightforward calculation we get the following system of equations for A' 



/ 




-1 


mi sin(2p) 


mi cos(2p) 





(1 — mi) sin(2p) 


(1- 


mi) cos(2p) ^ 




/ 


\ 









— 1— mi cos(2p) 


mi sin(2p) 





(mi —1) cos(2p) 


(1- 


-mi ) sin(2p) 













m 1 





-1 


1— mi 










• a = 




-iX 




(1- 


-7712) cos(2p) 


(1 — m 2 ) sin(2p) 





m2 cos(2p)— 1 


vri2 sin(2p) 














(1 


-ma) sin(2p) 


(m2 — 1) cos(2p) 





m2 sin(2p) 


— m2 cos(2p) — 1 














\ 










m-2 — 1 










/ 




V 


-iX J 



with coefficient vector a = (ai(l), a 2 (l), 03(1), ai(2), 02(2), 03(2)). Actually, there are two more equations for the 
variables ao(l) and ao(2), but these are already fixed to be zero in order to guarantee the condition m(g>~p (A 1 ) = 0. 
Solving these equations we get the result 



sm(p) 



1 - mi 



m 2 + mi 



1 — mi 2 — mi — m 2 



The comparison of sb(p) and sm(p) shows that the diffusive order of the quantum walk can distinguish between 
Bernoulli and Markov type decoherence, as opposed to the ballistic order which is the same for both control processes. 
In both cases, we see that in the coherent limit e — > and mi — > 1 the variance sb,m(j>) diverges independent of m-2. 
For the Markov process we can also consider the limit m 2 — > 1 which gives us localization, expressed by sm{p) —> 0. 

B. One dimensional Quantum Walks with equal position distribution 

In experimental implementations of one dimensional quantum walks one source of decoherence can be identified as 
dephasing of the internal degree of freedom of the walking particle [3J. Such a dephasing error can be modeled by 
introducing an additional z-Rotation R{9) of angle 9 before the application of the coin operation, such that a time 
step of the quantum walk Wg is then given by 



Weip) =C ■ R{6) ■ S = C 











s 



where C and S are the coin and shift operations as introduced in Eq. (25). In the following we will consider a 
quantum walk with Bernoulli type decoherence, the control process is assumed to chose a value 9 in each time step 
independently according to which Wg is applied subsequently. This quantum walk applied, up to special values of 9, 
Assumptions [T] [2] and [3] hence it is easy to compute the ballistic and diffusive scaling. 

One problem faced in experiments is that the quantum walk W = C ■ S can not easily be distinguished from its 
disturbed counterpart Wg = C ■ R(9) ■ S, because of the following little lemma. 

Lemma 10 For an initial state tp(p) — <j){z)e~ lzp localized at lattice point z the position distribution after any number 
of time steps t, generated by the quantum walks W = C ■ S and W = C ■ R{9) ■ S , are identical, where R(9) is an 
arbitrary z-Rotation. 



Proof In the Fourier picture it is easy to see that W and Wg deviate only by a constant momentum shift, i.e. 
Wg(p) = W(p + 9). Computing the state at lattice site x after t time steps, starting with a state ip initially localized 
at lattice site z we find: 



(W^Xx) = 7r ldp Wl{p)^p)e^ = — dp W t ( P + 9)^(z)e^-^ 



1 

2n 
1 

— < 
2tt 



I 



2?r 



ieix - z ^> / dp W t (p)t/j(z)e lp{x - z) = e'^-^iW'^ix) 
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Since both amplitudes deviate only by a global phase the probability to find the particle at a lattice point x, and 
therefore the position distribution, coincide. 

■ 

So if we start an experiment in a localized state we cannot distinguish between the two. On the other hand, however, 
implementing a quantum walk, that will in every time step apply the undisturbed walk W with some probability e 
and with probability (1 — e) apply the z- rotated quantum walk Wq, will for the most values of 9 result in diffusive 
behavior. 

As an explicit example consider as the undisturbed walk the one dimensional Hadamard walk as defined in Sect. 



IV A The two Kraus operators in momentum space are then given by 



*i(p) = HS = V 2 ( e iP _ e -iP ) KM = HR (°)S = V ^ ( -e~^ J (69) 

One can check that only the cases 9 g {0, tt, 2ir} lead to commuting Kraus operators that violate part three of 
Assumption [2] and are treated in example VIE In all other cases we get an asymptotic deterministic velocity v(p) = 
by Proposition [6j 




FIG. 7. (color online) The plot shows the variance s(p) for a dephased Hadamard Walk (691 with 9 = ^. 
Hadamard Walk is applied with probability e £ [0, 1] and the walk W-kIz is applied with probability (1 — e). 



The undisturbed 



To compute the diffusion constant we follow the instructions of Sect. VD and calculate the operator A' = ^ . a,j<jj, 
here decomposed in Pauli matrices. Choosing A! to have zero trace we find ao = and are left with three parameters 



a,j that can directly be inferred from Eq. ( 65 ) , which after eliminating equivalent conditions reads 



1 i(i + (i_e) e 2i(P+e) +ee 2i P) -((l- £ ) e 2 *(i>+ e > +ee 2i P) 
1 _i(i + (i_ e ) e -«(P+e> +ge -*j.) _((i- e ) e - 2 *(p+ e > +ee - 2 *p) 
10 -1 




Inferring A' from this equation, we can compute ~V'(A') and since 




(70) 



as in example VIA we can determine the variance from Eq. (61 ) 

A 2 cos(p + 6) 2 + e sin(2p + 6) sin(0) 



s(p) 



sin(<9) 2 



6(1 



(71) 



In contrast to example |VI A| we have an explicit p dependence of s. As one expects, in the limits e — 1 and e = 0, 
where only one coin is taken and the decoherence vanishes, the diffusion constant diverges. The same is true for the 
cases 9 € [0, tt, 2ir], where the Kraus operators commute and the quantum walk once again exhibits ballistic spreading 
as will be shown in example |VI E| So, also in the coherent limit 9 mod 2ir —> 0, where the two quantum walks, the 
Markov process chooses from, become equal, the diffusion constant diverges. 
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In Fig|7]we plot the variance s(p) for a z-Rotation of 6 = As explained in the last paragraph, we can observe 
divergence of s(p) for the coherent limits e = and e — 1. In between these two regimes the p dependence of s(p) is 
recognizable. 



C. Non-unitary Kraus operators 

The following example is a simple version of a quantum walk which satisfies Assumption[l]and[2]but not Assumption 
[3j This means there is no momentum transfer, and according to Proposition[5]the eigenvalue 1 of W is non-degenerate 
for almost all p, but the Kraus operators are non-unitary. Hence, the results of Sect. |V A| are applicable, and 
in particular formula (41) determines the asymptotic behavior of the expectation value of Q{t)/t, but now with a 
momentum dependent ve ocity v(p). 

We consider a one dimensional lattice with no internal degree of freedom. The decoherence will be of Bernoulli 
type, i.e. the Markov chain is actually trivial. The Kraus operators of the quantum walk are defined by 

{K^){x)= l -{^{x) + ^{x + l)) , {K 2 ^){x)= l -(^{x)-^(x-l)). (72) 
These operators satisfy the normalization condition K*Ki = 1 and their Fourier transforms are given by 

K 1 (p)= 1 -(l + e*n , K 2 (p) = ^(l-e- i n. 



In fact, the Kraus operators K\ and K% commute, hence, according to Sect. VIE the behavior of this quantum walk 
will be ballistic. The limit of the characteristic function of Q(t)/t is according to the discussion in Sect. Ill given by 

lim C Q[t)/t {\)= lim J dptr p(p)W[ /t (l) . 

Although our theory applies to this problem, it is instructive to calculate this limit using only Fourier methods. We 
need to compute the operator W^y t (l), in momentum space this becomes a one dimensional problem, for we have 

W £ (X) = J2 K* ( P )X(p)K l (p + Xe) = X(p) £ K* (p)tf,.(p + Xe) . 

i i 

But this means we can solve the equation W £ A £ — [i e A e exactly. Clearly, A £ = 1 and n e {p) = J2i ^* (p)Ki(p + Xe) 
solves the equation. Now, in ballistic scaling e = 1/t the operator W* / t (l) is just given by 

H* 1/t = j t ((1 + e- lp )(l +e ^ p+x ^) + (1 - e vp ){l - e -^ p+x/t) )) ' (73) 
' 1 + cos(A/t) + a(sin(p + X/t) - sin(p)) \ 1 ^ 

By a Taylor expansion of this equation to second order and using the formula lim t _ >00 (l + x/t + o(t -1 ))* = exp(x) 
once again we obtain the limit 

lim C Q{tyt (X)= [ d PP (p)e lXcos ^ 2 . 



As already pointed out, we would have obtained the same result by applying (41 ) and (17). And indeed, these formulas 
can be used to determine the asymptotic behavior of more general quantum walks. 

Proposition 11 Let the Kraus operators of a quantum walk in one lattice dimension with no internal degree of 
freedom be defined by 



{Kiip)(x) = a lk ip(x -k) O Ki(p) = Y a lk e ipk . 

k k 

where the coefficients have to satisfy the constraint X)i fc ^i(fc+a;) a ife — c o^o,x ^ x G 2 in order to guarantee the 
normalization condition ^2iK*Ki = 1. Then the group velocity and hence the asymptotic behavior of the quantum 
walk is determined by 

v(p) = Y \lx\ cos(px + arg^)) , j x = ^ o-i{k-x)a%kk ■ 



27 



Proof The statement follows directly from Eq. (41). 



Interestingly, we can infer from the characteristic function that the behavior of our example walk is truly ballistic. In 
fact, the second moment of Q(t)/t in the asymptotic limit is non-zero: 



A=0 



We would like to determine the asymptotic distribution for arbitrary initial states p and according to Proposition 
this requires knowledge of the group velocity v(p), which is for our example given by 



v(p) 



cos(p) 



By decomposing momentum space into subsets (— tt, 0] and (0, w] on which the function v(jp) is invertible and exploiting 
the point symmetry of v(p) we obtain the asymptotic distribution of Q(t)/t 



P(x) = 



y/1 - Ax 2 



where w _1 denotes the inverse of v restricted to (—tt, 0]. Again, we see that at the points where the derivative dv/dp 
vanishes, i.e. at the caustic points, there are peaks in the limiting distribution. If we choose p to be located at the 
origin, i.e. p(p) — 1/2tv, we get 



Ttyl — Ax 2 



The Fig. [8] compares this asymptotic distribution with the scaled probability distributions for a finite number of time 
steps. Our knowledge of the exact eigenvalue p e allows us to analyze the asymptotic distribution in further detail. 




The Taylor expansion of p £ to second order reads 

e 2 eX e 2 X 2 
p e = 1 +ep' + —p" + 0(e 3 ) = 1 + i— cos(p) — (isin(p) + 1) + 0(e 3 ) . 

Hence, the variance s(p) for this example is a complex valued function 

/ u2 // / U 2 2 + i2sin(p)-cos 2 (p) 2 
s(p)X = -p — (v - X) = A . 

The imaginary part of this "variance" defies its interpretation as the variance of added Gaussian noise. However, the 
corrections to the probability distribution computed from it will be real, after integration over momenta. In a similar 
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vein, ( 73 ) looks like the characteristic function of a sum of t random variables, but with complex "probabilities" . Only 



after integration over p the expression gives a probability distribution. 

Using the ideas from Section V C we can determine the 1 jt correction to the asymptotic distribution stemming 
from the second order of the perturbation expansion in e. Since we have A' = and we have chosen p(pi,P2) = l/2n 
we get 



a (A) 



i 

2^ 



dpe iXv{p) ( 1 - 



2 + i2 sin(p) — cos 2 (p) 
St 



A 2 ) +o(t- 1 ). 



With v{p) — cos(p)/2 the integral over p yields a sum of two Bessel functions for the first order approximation Ci(A, t) 
of C t (A), i.e. C* t (A) = d(A, t) + +o(^ 1 )' and 

<?i(V) = i ((8t - A 2 )J (A/2) - 2AJi(A/2)) . 

ot 



As already pointed out in Sect. IV B Ci is not integrable over R (cf. Fig. |9| and therefore we need to introduce a 



cutoff in the integration for the inverse Fourier transform of C t (A). In order to smoothen the resulting probability 
distribution we multiplied Ci(A,i) with a Gaussian g(X,t) = e~( x /') and computed the inverse Fourier transform 
numerically. This smoothening was necessary because of the rapidly oscillating behaviour of Ci(A,i), which can be 
seen in Fig. [9j leading to a poor convergence of the numerical integration. The resulting correction to the asymptotic 
position distribution for 10 time steps, also shown in Fig. [9j is in good agreement with the exact position distribution. 




D. Quantum walks with momentum shifts 

Now we want to drop Assumption [l] hence making the results of the preceding sections inapplicable. The quantum 
walk we are going to consider is closely related to the example of Sect. |VI C] the decoherence will still be of Bernoulli 
type and it has Kraus operators 

(K 1 <P)(p) = ±(l + e i ^)iP(p-q) , {K 2 iP)(p)= l -{l-e-^)iP{p-q). (75) 

The case q = recovers the quantum walk of the previous section which shows ballistic behavior. For g ^ we will 
use Fourier methods to show that the behavior is diffusive. For this purpose we will assume that q £ 2n ■ Q , i.e. 
there are two numbers m, n G Z relatively prime, such that q — 2ir ■ n/m. In position space these operators act in the 
following way 

(K^Kx) = l - (e ix «i>(x) + e ix «ip(x + 1)) , (K 2 ^)(x) = \ (e lx ^(x) - e ix ^{x - 1)) . 
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Although these operators themselves are not translation invariant because of the position dependent phase factor 
e lqx , the quantum walk constructed from them is translation invariant, as guaranteed by Corollary [3J And indeed, 
for arbitrary operators X the matrix elements of K*XKi pick up phase factors which depend only on the difference 
of the row and column index of the operator, which is clearly a translation invariant operation. 



Proposition 12 Let the Kraus operators be given by (75) with q € 2ir ■ Q. If q mod 2ir = the behavior is ballistic 



according to Sect. VIC For <7mod27r 7^ the ballistic order Q{t)/t converges to a point measure at the origin. In 
the case qmod2ir ^ it the diffusive order Q(t)/y/i converges to 

j/gmod27r = it the asymptotic distribution of Q(t) / ' \ft is given by 

P 2 (x) = ±J dXe-^e- 3 ^ J dpp(p)exp ( A 2 ^g^ 



Proof Similarly to the preceding example, we will calculate the quantity tr pW* (1) in the limit t — > 00 with appropriate 
scaling of e. We begin by noting that for an arbitrary function f(p,p') the following identity holds 



W e (f)(p,p') = ly^K^fe^Kje-^j (p,p') 

= 1(1 + e-*)(l + e^'+^)f(p + q,p' + q) + 
+ 1(1 - e lp )(l - e- l(p ' +£A) )/(P + Q,P' + ff) 

= ^ + cos (p' -p + e\)+i- (sm(p' + e\) - sm(p)))f(p + q,p' + q) . 

Exploiting the condition q = 2n ■ n/m, with n and m relatively prime, together with this identity we see that the 
characteristic function of the asymptotic distribution in e-scaling is given by 

/~ m—l ^ 

dpp(p)WUl)(p,p) — dpp(p) lim IT — (1 + cos (eA) + i • (sinCp+fcg+eA) - sin( P +fe 9 ))) s , (76) 
/ s— too A A 2 

J k=0 

with t = m-s. Hence, in ballistic scaling e — 1/t we get the result 



m— 1 

- A cos(p+kq) 



lim^C Q{t)/t (X) = / d P p{p) JJ e>™ 

■' fe=0 

dpp{p)e x ^=o ^ c ° s (P+ fc «) 

dppijp) 



1 



which follows from a Taylor expansion of (76) in e and lim t _ j . oc (l + a;/t+o(t ))* = exp(x). Therefore, all the moments 



vanish and the random variable Q(t)/t converges to a point measure at zero. On the other hand, some combinatorics 
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tells us that 

m—l 



JJ -(1 + cos(eA) + i ■ (sin(p + kq + eX) — sin(p + kq))) 

fe=0 

m_1 1 1 i 

[| ^(2 - ^e 2 \ 2 + i ■ cos(p + kq)e\ - - sin(p + kq)e 2 \ 2 + o(e 2 )) 

fe=0 

l_(2- + 2 ro - 1 feA5>os(p + fc< ? )-2 m - 3 £ 2 A 2 (^os(p + kq)\ + 

k \ k / 

+2 m-3 £ 2 A 2 ^ cos 2(p + fcg) _ m2 m-2 £ 2 A 2 _ j2"^ 2 e 2 A 2 ^ sill(p + kq) + o( £ 2 )) 
k k 

1 - m^e 2 X 2 + ^e 2 A 2 ^ cos 2 (p + kq) + o(e 2 )) 

k 

1 - ffe 2 A 2 + o(e 2 ) ,if gmod27r^7r 



1 - (3 - cos(2p)) +o(e 2 )) ,if gmod27r = 7r 



The first equality is just a Taylor expansion in s, in the second equality we sorted terms according to e and used the 
relation 

cos(p + kq) cos(p + Iq) = - f ^ cos(p + fcg) J - - ^ cos 2 (p + /eg) , 
k<l \ k J k 

and the third equality uses Y^T=o C0S iP + kq) — for m • gmod27r = 0. The last equality follows from the formula 



m— 1 m — l 



jr cos 2 (p + kq) = ^^2(l + cos(2 • (p + kq)) = ^ + <W cos(2p) . 

fe=0 fe=0 

Hence, for gmod27r ^ it and in diffusive scaling e = \j\ft we have for arbitrary states p 

f ( 3A 2 \ S f 3 * 2 3X '- 

lim C Q(t)/4 (A) = lim / dp p{p) 1 - — + o(.s" 1 ) = / dpp(p)e-^ = e~^s 

and the asymptotic probability distribution obtained via inverse Fourier transform reads 

_ , * 2 4a: 2 

P 1 (x) = —e-^. 

V07T 

If q mod 2-7T = 7r we get the result 

t limC Q(t)/t (A) = e-^ J dpp(p)ex P (\ 2C0S ^^ . 



In contrast to the case q = 0, the quantum walk for q ^ shows diffusive behavior. The asymptotic distribution 
together with the probability distribution for a finite number of time steps is shown in Fig. |10| 

E. Commuting Kraus operators 

In this subsection we want to give an example that violates the non-degeneracy condition of Assumption [2] and 
show that we can still compute the asymptotic position. We consider a decoherent quantum walk with Bernoulli type 
control process acting on a finite control space. The degeneracy of the eigenvalue 1 will be exactly the dimension of 
the coin space due to the assumption that all the Kraus operators Kj commute, but they must not have a common 
degenerate eigenspace. Since we will study uncorrelated noise in this example, we can forget about the classical control 
space of the Bernoulli process and model the decoherence by a quantum channel without explicit classical control. 
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ft *p(n/V7 ) 



ft *p(n/ft ) 



/ 



(a) 



(b) 



FIG. 10. (color online) Asymptotic position distribution for a walk (75 I with momentum shift q = 7r/16 (green/dotted line) 



after (a) 20 steps and (b) 100 steps. The red/dashed line is the asymptotic distribution 



Proposition 13 Suppose that a quantum walk W with Bernoulli type decoherence on a finite control space is im- 
plemented with normal and commuting Kraus operators Kj(p) without degenerate eigenvalues. The Kraus operators 
can then be decomposed in a common eigenbasis Kj(p) = kj ta (p)\^ a (p)}(^ a (p)\. IfJ2j kj^ a kj,p ^ 1 for a ^ (3 the 
asymptotic position distribution in ballistic scaling is given by 



with the vector operator 



t lim C Q ( t )/ t (A) = dp tr(p(p)< 



i\-V 



V T = kj,g{p) 



dp T 



\Mp))(Mp)\ ■ 



Proof Since all the Kraus operators K 1 (p) commute and are normal by assumption, they can certainly be diagonalized 
in a common eigenbasis 

K j(p) =^ k 3,<x(p)\i>ot(p)){i><x(p)\ ■ 



One can check that a basis of invariant operators of W is given by the one-dimensional projections \ip a (p)){ij)p(p)\. 
Since W is a quantum channel it follows that \kj, a \ 2 = 1 and together with the condition on the sums of eigenvalues 
for a ^ j3 this implies that W(j4) = A holds only for operators in the span of {\ip a (p)) (ip a (p)\ — P a }. 

In order to compute the ballistic scaling we arc interested in the corrections to the eigenvalue one. Evaluation of 
the first perturbation order on an arbitrary operator R — ^ Q r a P a in the span of the eigenoperators to the eigenvalue 
one yields 

W(R'){p) — R' = n'R - W'(R) 
If evaluated with respect to the state pp = \ipp(p))(i/jp(p)\ the left hand side of this expression vanishes and we have 



'V 



Q.J 



£ = 



Y k iAp)^ k iAp) ■ A + Y k i,p(p) k J,<*(p) tr [pp ~d~ £ Pa ( p + eA ) 



J. a 



e=Q 



where we just used the definition of R and K'^ip) and the fact that {IV'a)} is an eigenbasis of the Kraus operators. 
Since Pp(d p P a )Pp — holds for orthogonal projectors, the second summand on the right-hand side vanishes and we 
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get 

3 

Substituting these results into W 1( / f (l) we find 

fim W* /t (I)(p) = lim^(l + ^ a {p))*P a =J2 e "' a(P)p «(P) 



Together with ( 10 1 and ( 12 ) we have 

Urn C Q(t) (X) = J2 J d P tr (p(p)^ ip) P a ( P )) 



which ends the proof. 



Corollary 14 For a unitary implemented walk, fulfilling the conditions of Proposition 13 the components of the 
operator V are given as the weighted sums of the components of the group velocity operators Vj of the single walk 
operators Wj 



where r\j is the probability that the quantum walk Wj is applied in a time step. 

Proof The unitarity of all the Kraus operators Kj implies that the eigenvalues fcj, a (p) are given by phases ^/fjJe iuJ: '' a ^ p ' , 
r)j being the probability to apply walk operator Kj. This implies kj^ a (p) = Vjkja(p), and therefore we get 

dkj ia (p) . du j>a 



k 3,a(P) ^ = ir lj 



dp T 8p T 

Inserting this result into the definition of V T finishes the proof. ■ 

To conclude this subsection we will look at a one dimensional example. For a given unitary one dimensional 
quantum walk W(p) we chose the two Kraus operators 

K 1 (p)= l -{1 + W{p)) K 2 (p) = 1(1 - W(p)) . 

Since both Kj are just functions of the original walk operator W(p) they will certainly commute, and since W is 
unitary and therefore diagonalizable, so are the Kj with eigenvalues 

K±{p) = \{t + ^ ±(p) ) *2,±(p) = \i l - ^ ±{p) ) ■ 



One can calculate that for oj+(p) 7^ u_{p) the kj t ± satisfy the sum-condition of Proposition 13 
compute V directly from the eigenvalues 

j.a a 



and therefore we can 



where P a are the eigenprojections and Vw the group velocity operator of the original quantum walk W . So this 

dp 



decoherence model halves the velocities ^ of the undisturbed quantum walk. 
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